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We numerically investigate the Olami-Feder-Christensen model for earthquakes in order to char- 
acterise its scaling behaviour. We show that ordinary finite size scaling in the model is violated 
due to global, system wide events. Nevertheless we find that subsystems of linear dimension small 
compared to the overall system size obey finite (subsystem) size scaling, with universal critical coef- 
ficients, for the earthquake events localised within the subsystem. We provide evidence, moreover, 
that large earthquakes responsible for breaking finite size scaling are initiated predominantly near 
the boundary. 



PACS numbers: 05.65. +b, 45.70.Ht 
I. INTRODUCTION 

Many dynamical phenomena in nature are intermit- 
tent. This "bursty" dynamics may be related to an un- 
derlying complex state, often characterised by long range 
correlations in space and time. For example, the crust of 
the earth alternates long periods of relative quiescence 
with burst of activity (earthquakes), which have a wide 
range of possible sizes. The behaviour of earthquakes 
is described by the empirical Gutenberg-Richter (GR) 
law 101, where the distribution of energy dissipated in 
earthquake events is a power law over many orders of 
magnitude in energy. The GR scaling extends from the 
smallest measurable earthquakes, which are equivalent to 
a truck passing by, to the most disastrous that have been 
recorded. Similar scale free behaviour of bursts is ob- 
served in vastly different kinds of physical systems such 
as flux motion through disordered type II superconduc- 
tors placed in a mag netic field §, or in granular piles, 
under some conditions etc. 

Self-organised criticality (SOC) Q has been proposed 
as a general dynamical principal behind the observed 
complex behaviour of many natural phenomena. It refers 
to the fundamental property of slowly driven, extended 
systems to organise, over a sufficiently long transient pe- 
riod, into a dynamical critical state which lacks any char- 
acteristic time or length scale. The amplitude of the re- 
sponse of the system to an external perturbation follows 
a power law distribution. A number of simple models 
have been developed to test the appHcabihty of SOC 
to a variety of complex interacting dynamical systems, 
such as sand piles and earthquakes (for a review see e.g. 
ref. 

One of the basic theoretical problems is to identify ro- 
bust, and thus physically relevant mechanisms for SOC to 
emerge, and to define a space of parameters and dynam- 
ical processes where SOC is a stable feature. Although it 
has been proposed that the presence of conservation laws 
(e.g. sand grains being transported in a sand pile) or spe- 
cial symmetries was necessary for SOC |^,||, many known 
examples of physical phenomena and some models have 



been found where no apparent conservation law or special 
symmetry exists. These include, besides earthquakes, bi- 
ological evolution, forest fires, epidemics, (possibly) solar 
flares, (possibly) reconnection events in the magnetotail, 
etc. ||^-|^. In contrast to conservative systems, the mech- 
anisms for SOC in non-conservative systems are not very 
well established. 

A nonconservative SOC model that in recent years has 
attracted much attention is the so called OFC model . 
The OFC model is a simplified lattice representation of 
a spring-block model for earthquake dynamics which was 
originally introduced by Burridge and Knopoff pd| ] . The 
Burridge-Knopoff model can be schematized as a two di- 
mensional network of blocks interconnected by springs. 
All blocks are subject to an external driving force, which 
pulls them, and to a static friction, which opposes their 
motion. In the OFC model, each site of a lattice is as- 
sociated with a continuous variable, which represent the 
force acting on a block. A slow driving is applied to 
the system by increasing uniformly and simultaneously 
the forces of all the elements. When the force on a site 
exceeds some threshold value (the maximal static fric- 
tion), the site relaxes and distributes part of its force to 
nearest neighbour sites. Each such discharge event is ac- 
companied by a local loss in accumulated force from the 
system. This conceptually simple and seemingly numeri- 
cally tractable model reproduces some of the qualitative 
phenomenology of the statistics of earthquake events such 
as power law behaviour over a range of sizes, and inter- 
mittency or clustering of large events iQ . Extensions of 
the model have been recently developed which reproduce 
to some extent Omori's law and other temporal patterns 
associated with earthquakes . 

In the context of nonconservative models, the OFC 
model is of particular interest as it is possible to directly 
control the level of conservation of the dynamics through 
a parameter a. Early analysis on relatively small sys- 
tems indicated that the OFC model exhibited SOC, in 
the sense that earthquakes in the steady state obeyed 
finite size scaling (FSS) when the system size was var- 
ied pO|. However, the critical coefficients obtained using 
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the FSS ansatz were found to be nonuniversal. In par- 
ticular the exponents characterising the power law dis- 
tributions appeared to vary with both the dissipation 
parameter, a, and the form of the boundary conditions. 
This would have been in sharp contrast to the usual fixed 
point picture of critical phenomena where most micro- 
scopic details are irrelevant. Moreover some apparent 
critical exponents obtained using FSS violated physical 
bounds , putting some doubt on the existence of criti- 
cality in the model. Recently it was shown using a multi- 
scaling analysis of large-scale simulations that, actually, 
the avalanche size distribution has a universal power law 
behaviour, independent of the dissipation parameter and 
for different boundary conditions, but that the cutoff in 
the power-law distribution does not behave according to 
FSS ll^. In larger systems, proportionally more of the 
force can be dissipated in the largest events that occur, 
and the cutoff function becomes sharper and sharper as 
the system size increases. 

Departures from standard FSS have been reported 
for various SOC models as, for example, the sand- 
pile modeljlfl], the Drossel-Schwabl (DS) forest fire 
model and the Zhang model In this paper we 

address the question of the origin of the breaking of FSS 
and its relation to the mechanism responsible for SOC 
in the OFC model. In particular we test the implicit 
assumption behind the FSS hypothesis that a finite sys- 
tems behaves as a subsystem of a larger system. The 
paper is organised as follows. In the section II we de- 
scribe in some detail the model. In section III we present 
the results of our numerical study relative to two differ- 
ent type of probability distributions for earthquake sizes. 
The first distribution concerns earthquakes which are lo- 
calised within a given subsystem. We show that this 
subset of earthquakes exhibits ordinary FSS as long as 
the linear extent. A, of the subsystem is sufficiently small 
compared to the linear extent, L, of the entire system. 
The second distribution groups earthquakes according to 
the position of their starting site relative to the bound- 
aries of the system. From this investigation, we deduce 
that FSS is violated due to large events initiated in a re- 
gion near the boundary. Finally, in section IV we discuss 
our results and draw some conclusions. 



1. Uniform drive: all forces Fi are increased at the 
same rate, until one of them reaches the value Fth- 

2. Earthquake: when a site becomes unstable (i.e 
Fi > Fth), the uniform driving is stopped and the 
system evolves according to the following local re- 
laxation rule 

F.>Fth^\''^^^ (1) 

until there are no more unstable sites. In eq. (|^), 
the subscript "nn" stands for the four nearest 
neighbours to site i. 

Since only a fraction, Aa, of the force is redistributed in 
each relaxation event (toppling) , the model is nonconser- 
vative for a < 1/4. In the following we concentrate on 
this case, i.e. < a < 1/4. 

To completely define the model we need to specify the 
boundary conditions. Boundaries are believed to play a 
crucial role for the observation of critical behaviour in 
the OFC model. It has been suggested that they act as 
inhomogeneities which frustrate the natural tendency of 
the model to order into a periodic state ||20| , ^ . Indeed, 
for sufficiently small values of the conservation parame- 
ter a (a < ac — 0.18), a system with periodic boundary 
conditions quickly reaches an exactly periodic state with 
only earthquakes of size one. For larger values of a the 
situation is slightly more complicated with multiple top- 
plings involved in a single avalanche, but the avalanches 
are still localised and criticality is not observed ||2^. A 
system with open boundaries is prevented from reach- 
ing a periodic state because boundary sites have fewer 
neighbours and therefore cycle at a different frequency 
from bulk sites. Middleton and Tang suggested that 
the inhomogeneity created by the boundaries propagates 
into the bulk of the system, developing, in this way, long 
range spatial correlations. They named this mechanism 
"marginal synchronisation" or phase locking [ p2[ . In ac- 
cordance with previous studies, therefore, we consider 
open boundary conditions. If a boundary (corner) site 
topples, an extra amount aFi {2aFi) is simply lost by 
the system. 



II. THE MODEL 

We consider a two-dimensional square lattice of L x L 
sites. To each site i of the lattice we associate a continu- 
ous variable Fi , which initially take some random values 
between zero and a threshold value Fth- The dynamics 
proceeds then indefinitely. In the limit of infinite time 
scale separation between the slow driving and the (al- 
most) instantaneous earthquake process, the dynamics 
is: 



III. RESULTS: PROBABILITY DISTRIBUTIONS 

After a sufficiently long transient time, the system 
reaches a stationary state. Several statistical properties 
can be used to characterise this state. Most previous 
studies of the OFC model have focused on the behaviour 
of the probability distribution of earthquake sizes, Pl(s), 
where L is the size of the system and s is the total number 
of topplings events during an earthquake p0| , p^ , p0| -p5[ . 
We choose instead to analyse the behaviour of different 
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distributions for avalanches sizes, which distinguish be- 
tween earthquakes according to the region of the lattice 
involved (e.g. bulk or boundary) and the coordinates 
of the triggering site (see figure |]). This investigation 
is particularly relevant for the OFC model in view of 
the strong inhomogcncity in the spatial distribution of 
avalanches [ p0| , p2[p5| ]. According to ref. |^^, for exam- 
ple, large avalanches are localised near the boundary (at 
least for a < Q!c). As a minor technical remark, we point 
out that we exclude from our data avalanches which in- 
volve only one site (s = 1) as they appear to behave 
according to their own statistics [ pO| . As we are mainly 
interested in asymptotically large earthquakes, this does 
not alter our conclusions. 

We consider a subsystem of linear extent A centred 
in a system of size L. The first distribution we in- 
troduce, Pconf{.^,L,s), is the normalised distribution of 
earthquake sizes restricted to earthquakes which are con- 
fined entirely within the subsystem (e.g. avalanche (a) in 
fig. 0). The model is driven according to its usual dynam- 
ics but only those particular earthquakes are counted. 
According to our definition, the case X = L corresponds 
to the distribution of avalanches which do not reach the 
boundary of the system. As shown in fig. ^, the dis- 
tribution PconfiX,L,s) becomes independent of L, if L 
is considerably larger then A (approximately L > 2A). 
When L approaches A, this is no longer the case and the 
cutoff in the distribution is pushed to larger sizes. Al- 
though we have shown in Fig. ^ only the distributions for 
a = 0.18 and A = 32, 64, analogous considerations apply 
to different values of a and for different sizes, A. Since 
for a generic L, Pconf{.\L,s) ^ PconfiX,\s), a small 
portion of a large system is substantially different from a 
finite system of the same size, contrary to what happens 
in equilibrium critical phenomena. A similar observa- 
tion was made in ref. Ijl^l for the DS forest fire model. 
In the following, we denote with Peon/ (A, s) the distri- 
bution Pconf{\ L, s) in the limit where the distribution 
does not appear to depend on L. In order to determine 
numerically these distributions, for each value of A we 
have simulated (for at least 2 • 10^ earthquakes) a system 
of size L = 2A. The dependence on L of Pconf can in this 
case be safely neglected. With this choice, the accuracy 
of the measures and the range of scales investigated are 
optimised, within our computational limits. 

In Fig. II we report a FSS collapse of Pconf {X, s) for dif- 
ferent values of a . Contrary to the entire distribution of 
earthquake sizes, Pl{s), we observe that Peon/ (A, s) sat- 
isfies the FSS hypothesis, i.e. Peon/(A, s) ~ f{s/ X^), 
with universal critical coefficients. The curve correspond- 
ing to a = 0.15 and A = 256 shows some noisy behaviour, 
due to the difhculties in collecting good statistics in this 
case. Indeed by decreasing a, the relative fraction of 
earthquakes in the bulk of the system (with size s > 1) 
diminishes. Nonetheless, there is no evident sign that 
FSS is violated in this case. The critical exponents used 



in the fit of Fig. ^ are /? = 3.6 and D = 2, independent of 
the dissipation parameter a. The value of the histogram 
exponent t — P / D 1.8 we obtain is the same as that 
found for Pl{s) ||15|. In addition, the value of D we find 
corresponds to the largest possible value for the entire 
distribution {Dmax in ref. [p5[), as it can be shown that 
non-conservation requires D < 2 . 

The scaling behaviour of Pconf appear to be reasonably 
robust with respect to translating the subsystem within 
the entire system; in fig. ^ we report a FSS plot for the 
subsystem placed on a boundary and on a corner of the 
system for a — 0.18. While the FSS collapse for the sub- 
system placed on the boundary is rather good, some devi- 
ations from FSS are observed in the cut-off region for the 
case of the subsystem placed in the corner. We believe 
this behaviour can be ascribed to the enhanced bound- 
ary effects in the latter case (two side of the subsystem 
are boundary sides instead of only one) and would dis- 
appear if larger (sub)systems could be studied. This pic- 
ture is confirmed by choosing different a values: for the 
subsystem in the corner, deviations from FSS are more 
pronounced for a = 0.21 and are absent for a = 0.15. 
For the subsystem on the boundary, instead, the quality 
of FSS collapse is rather convincing in all cases. 

We introduce next the distributions P< (A, L, s) and 
P>(A, P, s). These are the normalised distribution of 
earthquakes which start respectively within (P<) and 
outside (P>) the subsystem of size A, irrespective of 
whether they stay in or go out of the subsystem (see 
fig. 1^). The only difference between these two distribu- 
tions, P< and P>, is the location of the site that triggers 
the avalanche. We observe numerically that the distri- 
butions P<(A, P,s) and P>(A, P, s) become independent 
of A respectively in the limit A ^ P and A ~ P. As an 
example we report in Fig. |5|the behaviour of P<(A, P, s) 
and P>(A, P, s) for a = 0.18, L — 256 and for various A. 
For simplicity, in the following we denote with P^{L, s) 
and P>{L, s) the distributions in the limit where they do 
not depend on A. 

We consider therefore two centred subsystems of linear 
extent A2 > Ai , such that the above conditions are satis- 
fied. More specifically, we choose Ai = ^P and A2 = |P. 
In this case P> corresponds to the subset of earthquakes 
which are initiated in some "boundary" region and P< 
corresponds to the subset of earthquakes which are ini- 
tiated within some "bulk" region. In fig. ^ we report a 
FSS scaling plot both for P> and P<. In this figure we 
choose P = 2 as the maximum allowed value. It is clear 
that the boundary distribution P> cannot be collapsed 
according to the FSS ansatz. In fact it develops a sharper 
and sharper cutoff which changes shape and which has 
an excess of large events (the cutoff moves towards right 
for increasing P). The bulk distribution P< instead does 
not develop a noticeably sharper cutoff and does not ap- 
pear to change its shape. It may possibly be collapsed 
according to the FSS ansatz. Consistent with the results 



3 



for Pconf and for P^, the power law exponent for P< is 
T = /3i/D ~ 1.8. 

The above numerical analysis indicates that the large 
events which violate FSS are triggered by sites in a 
boundary region. Indeed the behaviour of the cutoff for 
the collapsed probability distributions P> and Pl is very 
similar (see Fig. 1 in ref. [^). Although large earth- 
quakes are focussed mainly toward the boundary, as sug- 
gested in ref. ||2^, they occur also in the bulk of the 
system, even for low values of a, as can be deduced from 
Fig. ^ and Fig. ^. Moreover, we do not observe any sig- 
nificant qualitative change in the behaviour of the system 
around a = ttc — 0.18, as claimed in ref. [ pO| . 

IV. DISCUSSION AND CONCLUSIONS 

Similarly to other SOC systems ||l^[l9|], the nonconser- 
vative OFC model shows relevant deviations from simple 
FSS ll^. In this paper we have investigated the origin 
of this phenomenon, finding that FSS in the OFC model 
is violated because of large, system wide earthquakes. In 
fact, we have found that earthquakes localised within a 
given subsystem do obey ordinary FSS, with universal 
critical exponents, independently of whether the subsys- 
tem is placed at the centre or on a boundary of the sys- 
tem. The value of the power law exponent, r ~ 1.8, 
for the "confined" distribution agrees with the one for 
the entire distribution. We have shown, moreover, that 
the probability distribution for earthquakes initiated in 
a boundary region do not obey FSS, because of an "ex- 
cess" of large events. This would result in an apparent 
exponent D > 2 which is not allowed in the nonconser- 
vative case. On the other hand, the probability distribu- 
tion for earthquakes starting in the bulk of the system 
is compatible with a FSS hypothesis. In particular, the 
critical exponent is Z? = 2 in this case, indicating that 
large earthquakes responsible for breaking finite size scal- 
ing are initiated predominantly near the boundary of the 
system. 

Self-organised criticality in the OFC model has been 
ascribed to a mechanism of "marginal synchronisa- 
tion" A system with open boundaries becomes al- 
most synchronised by an invasion process where spatial 
correlations develop from the boundaries. It was sug- 
gested that sites close to the boundaries start to organise 
themselves first, building up long range correlations. The 
critical region grows with time, until, in the stationary 
state, it invades the whole lattice . Our findings on the 
large events occurring at the boundary seem to indicate 
that the effect of synchronisation is stronger for boundary 
sites than for bulk sites. This view is supported also by 
the "on screen" observation that large earthquakes tend 
to be triggered repetitively by the same sites over a long 
time scale (a result which seems to be confirmed also by 
the study in ref. p5|). 
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L, the dashed hne the subsystem of Mnear extent A. Triggering ^ | ^ ^ ^ 

sites arc denoted with a full circle, toppling sites with a cross. -6 -4 -2 

Avalanche (a) contributes to the distribution Pcon/(A, L, s); log,„(s/X°) 
(a) and (b) to P<(A,L,s); (c) to P>(A,L,s). 

FIG. 3. Finite-size scaling plots of Pcon/(A,s) (with the 
subsystem placed at the centre) for (a) a = 0.15, (b) a = 0.18 
and (c) a = 0.21. The critical exponents are /3 = 3.6 and 




-6 -4 -2 



FIG. 4. Finite-size scaling plots of Pconf for the subsystem 

placed (a) on a corner and (b) on a boundary of the system 
{a — 0.18). The critical exponents are /3 = 3.6 and D = 2; 
the slope of the straight line is r = 1.8. Curve (a) has been 
shifted, X —f X — 1. 
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FIG. 5. Probability distributions (a) P<(A,L,s) and (b) 
P> (A, L, s) for a = 0.18, for L = 256 and for various A. 



10 




_g I \ , \ , \ , L_ 

-6 -4 -2 

iog,o(s/L°) 



FIG. 6. Finite-size scaling plots of (a) P< (L, s) and (b) 
P>(L, s). Different sets of curves corresponds, from left to 
right, to a = 0.15 (shifted by a; a; — 1), a — 0.18 and 
a = 0.21 (shifted by a; — > 1). The exponents are /3i = 3.6, 
P2 = 3.9 and D = 2. 



6 



